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Abstract 

We  consider  an  optimization  problem  for  a  multiclass  queueing  model  of  a  de- 
centralized production-distribution  system:  A  single  server  can  produce  a  variety  of 
products,  and  completed  units  pass  through  an  infinite  server  delay  node  on  their 
way  to  a  warehouse  inventory.  Each  product  has  its  own  warehouse,  demand  for  each 
product  is  satisfied  by  its  associated  inventory,  and  unsatisfied  demand  is  backordered. 
Each  warehouse  operates  independently  under  its  own  base  stock  policy:  The  initial 
inventory  is  set  at  a  base  stock  level,  and  each  unit  of  demand  triggers  an  order  with 
the  server.  The  server  processes  the  orders  in  a  first-come  first-served  manner,  and 
the  optimization  problem  is  to  find  the  base  stock  level  for  each  product  to  minimize 
the  expected  average  inventory  holding  and  backordering  cost.  A  heavy  traffic  approx- 
imation of  this  problem  yields  an  inventory  process  that  possesses  both  diffusion  and 
Gaussian  components.  Computational  results  are  performed  to  cissess  the  accuracy  of 
our  analysis. 


1      Introduction 

The  travel  time  between  nodes,  or  queues,  in  a  queueing  network  model  is  typically  assumed 
to  be  instantaneous.  However,  customers  in  most  real  queueing  systems  incur  delays  when 
traveling  from  one  node  to  another.  When  these  delays  are  essentially  independent  of  the 
congestion  in  the  network,  as  is  often  the  case,  the  delays  are  appropriately  modeled  as 
infinite  server  nodes  within  the  queueing  network.  Examples  include  the  delay  for  a  vehicle 
to  travel  between  traffic  intersections,  a  job  to  move  between  workstations  in  a  factory,  and 
data  to  be  transmitted  over  a  dedicated  low  speed  transmission  line.  An  infinite  server  node 
may  also  be  used  as  an  aggregate  representation  of  a  large  number  of  non-bottleneck  (that 
is,  low  utilization)  workcenters  in  a  large  factory,  such  as  a  wafer  fabrication  facility.  The 
heavy  traffic  approach  developed  by  Harrison  (1988)  has  been  used  in  recent  years  (see, 
for  example,  Wein  1992a,  Kelly  and  Laws  1993  and  references  therein)  to  control  queueing 
networks  that  possess  no  infinite  server  nodes.  The  aim  of  this  paper  is  to  make  a  first  step 
towards  generalizing  this  approach  to  queueing  networks  that  contain  infinite  server  nodes. 

Towards  this  end,  we  focus  on  a  specific  problem  that  is  simple  but  well  motivated:  An 
idealized  queueing  model  of  a  multi-product  production-distribution  system  that  consists  of 
a  single  production  facility  and  K  warehouses.  The  facihty  produces  K  different  products, 
and  a  completed  unit  of  a  given  product  is  transported  to  its  respective  warehouse;  hence, 
all  inventories  are  held  at  the  warehouses.  These  inventories  service  the  actual  customer 
demand,  and  unsatisfied  demand  is  backordered.  In  our  queueing  model,  the  production 
facility  is  modeled  as  a  single  server,  and  each  product  has  its  own  general  service  time 
distribution  and  renewal  demand  process.  The  delay  incurred  by  a  unit  to  travel  from  the 
facility  to  a  warehouse  is  a  random  variable.  As  alluded  to  above,  one  could  also  interpret 
the  model  as  a  production-inventory  system  with  a  single  bottleneck  workstation,  where  the 
delay  represents  the  aggregate  sojourn  time  through  the  workstations  downstream  of  the 
bottleneck. 

The  performance  of  such  a  system  relies  heavily  upon  the  production  control  policy  that 
is  employed.  We  derive  a  control  policy  for  a  decentralized  system:  Each  warehouse  has  an 


initial  base  stock  level  of  inventory,  and  whenever  a  demand  for  a  unit  of  product  i  occurs,  the 
inventory  at  warehouse  i  is  depleted  and  the  warehouse  places  an  order  with  the  production 
facihty.  The  facility  serves  orders  in  a  first-come  first-served  (FCFS)  fashion  and  is  idle  if 
no  orders  are  in  its  queue.  The  optimization  problem  for  the  decentrahzed  problem  is  to 
derive  a  base  stock  level  for  each  warehouse  to  minimize  the  average  cost  of  holding  and 
backordering  inventory. 

Since  the  optimization  problem  appears  to  be  difficult  to  analyze  without  making  addi- 
tional distributional  assumptions,  we  resort  to  a  heavy  traffic  analysis.  This  approprimation 
procedure  requires  the  server  to  be  busy  the  great  majority  of  time  to  meet  average  demand, 
and  assumes  that  the  travel  delays  are  lengthy.  Under  these  assumptions,  the  inventory  pro- 
cess for  each  product,  appropriately  normalized,  has  a  diffusion  and  a  Gaussian  component. 
Although  we  have  been  unable  to  derive  the  stationary  distribution  of  this  inventory  pro- 
cess, our  analysis  reveals  the  structure  of  the  hmiting  behavior  of  the  production-distribution 
system.  Under  additional  assumptions,  we  employ  our  analysis  to  numerically  derive  base 
stock  levels,  and  find  that  they  are  very  close  in  value  to  the  optimal  levels  obtained  via 
simulation. 

We  initially  planned  to  also  analyze  a  centralized  system,  where  a  controller  has  global 
information  about  the  inventories  of  each  product,  and  decides  in  a  dynamic  fashion  which 
product,  if  any,  to  produce  next.  Comparison  of  these  two  systems  would  allow  us  to  assess 
the  value  of  centraHzed  information  in  production-distribution  systems;  see  Zheng  and  Zipkin 
(1990)  and  Wein  (1992b)  for  similar  comparisons  for  production-inventory  systems,  which 
are  essentially  production-distribution  systems  with  no  travel  delays.  During  the  course  of 
our  analysis,  we  developed  a  variant  of  a  myopic  policy  proposed  by  Pena  and  Zipkin  (1992). 
Although  numerical  results  in  Pena  and  Zipkin,  Veatch  and  Wein  (1993)  and  Krichagina  et 
al.  (1992)  show  that  this  policy  performs  well  in  a  production-inventory  setting,  we  found 
that  our  variant  of  the  policy  did  not  perform  well  for  the  production-distribution  system; 
in  fact,  the  system  did  not  stabilize  under  this  policy  in  our  simulation  study.  It  appears 
that  the  long  travel  delays  lead  to  a  system  that  is  not  easily  controllable.    Consequently, 


this  paper  does  not  contain  an  analysis  of  the  centralized  problem. 

In  summary,  control  policies  for  queueing  networks  with  infinite  server  delay  nodes  are 
very  difficult  to  analyze;  if  the  delays  are  not  exponential,  then  the  cardinality  of  the  state 
space  can  become  unwieldy.  Unfortunately,  our  analysis  suggests  that  these  types  of  problems 
will  not  yield  easily  to  heavy  traffic  analysis. 

The  remainder  of  this  paper  is  organized  as  follows.  The  queueing  control  problem  for 
the  decentralized  model  is  defined  in  Section  2,  and  a  sequence  of  systems  approaching 
heavy  traffic  is  introduced  in  Section  3.  An  approximating  fluid  control  problem  and  an 
approximating  diffusion-Gaussian  control  problem  are  formulated  in  Section  4;  the  fluid 
analysis  is  a  prerequisite  to  the  more  elaborate  heavy  traffic  approximation,  and  allows 
the  controls  in  the  heavy  traffic  model  to  be  expressed  as  deviations  from  the  optimal  fluid 
controls.  Further  analysis  of  the  model  is  undertaken  in  Section  5,  along  with  a  computational 
study  that  assesses  the  accuracy  of  our  results. 

2      The  Problem  Formulation 

The  decentralized  production-distribution  system  pictured  in  Figure  1  consists  of  a  single 
machine  producing  K  products.  Units  of  product  i  =  I, . . .  ,K  have  a  general  service  time 
distribution  with  service  rate  //,  and  squared  coefficient  of  variation  i/^,.  We  let  Ai{t)  = 
max{j  :  ^1  +  •  •  ■  +  ^j  <  0  denote  the  renewal  process  associated  with  the  iid  service  times 
{i^j}i>i-  Completed  units  incur  a  random  delay  before  entering  the  finished  good  inventory 
at  the  appropriate  warehouse.  Let  {i]j}j>i  be  the  sequence  of  iid  delays  for  product  i,  and 
define  the  continuous  distribution  function  Ff{t)  =  P{ti[  <  t).  We  denote  the  mean  of  this 
distribution  by  A~^.  Each  product  has  an  independent  renewal  process  Di{t)  =  max{j  : 
C\  +  ■  •  ■  +  Cj  ^  0'  which  is  the  number  of  product  i  units  demanded  up  to  time  t.  We  let 
di  denote  the  demand  rate  and  i/J,  be  the  square  coefficient  of  variation  of  the  interdemand 
times. 

We  assume  that  warehouses  operate  independently.  Each  warehouse  initially  has  C, 
units  of  inventory,  and  each  unit  of  demand  for  product  i  simultaneously  triggers  two  events: 


Product  Vs  inventory  is  depleted  and  a  production  order  is  placed  with  the  machine.  The 
machine  serves  orders  in  a  FCFS  fashion  and  is  idle  if  no  orders  are  in  the  queue.  If  we  let 
Tj{t)  represent  the  total  lime  in  [0,/]  devoted  to  the  production  of  product  z,  then  .4  i/^/)) 
is  the  total  numbe'-  of  product  i  units  produced  up  to  time  t.  Denote  l>y  Q,{t)  the  number 
of  product  i  orders  in  queue  or  in  service  at  the  machine  at  time  t.  Let  Zt{t)  be  the  r.r,'.Qber 
of  unite-  in  product  i's  delay  node.  Then 


Q.(0  =  A(0  -  MT,{t)) 


[•2.r 


and 


(2.2) 


where  r'  is  the  time  of  the  j'^  jump  of  the  process  A,{T,{-));  it  corresponds  to  the  epoch 
when  the  j*-^  unit  of  product  i  completes  production. 
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Figure  1.  The  production-distribution  system. 

Let  A'.(f)  den.jte  the  number  of  units  in  product  i's  inventory  at  time  \.  This  quantity 
will  be  negative  i'  units  aie  backordered,  and  the  vector  A'(-)  =  {X^{-))  will  be  called  the 
inventory  process,  'jiven  -V,(0)  =  Cj,  the  inventory  process  can  be  expressed  as 


X.(0  =  c,  -  QAt)  -  z^t)  ,   t>0  . 


(2.3) 


Let  hi  and  6,  be  the  respective  costs  per  unit  time  of  holding  and  backordering  a  unit  of 


product  i,  and  define  the  cost  function 

H{x)  =  ^H,ix), 
1=1 

where  Hi{x)  =  hixf  +  biX~ ,  x'^  =  max(x,0)  and  x~  —  min(a:,0).  The  production  control 
problem  is  to  choose  Ci , . . . ,  Ck  to  minimize  the  expected  average  cost 

J(Ci, . . . ,  Ck)  =  lim sup  ^E  f  H{X{s))ds  .  (2.4) 

T— oo      J         ^0 

Since  the  optimization  problem  (2.1)-(2.4)  does  not  appear  to  be  tractable  without  additional 
distributional  assumptions,  we  turn  to  the  aid  of  asymptotic  analysis  and  assume  that  the 
system  operates  in  heavy  traffic;  that  is,  the  traffic  intensity 

P  =  E^<1  (2.5) 

«=i  f^' 

is  close  to  1.    Assumption  (2.5)  requires  that  the  machine  be  working  the  great  majority 

of  time  to  satisfy  demand  over  the  long  run.  Hereafter,  the  production-distribution  system 

described  in  this  section  will  be  referred  to  as  the  original  system  S. 

For  subsequent  use,  we  introduce  the  large  parameter 

no  =  (1  -  />)-'  (2.6) 

and  the  rescaled  distribution  functions  of  delay  times 

F.{t)  =  Ff{t{l-p)-')^F^{not).  (2.7) 

3      A  Sequence  of  Systems  in  Heavy  Traffic 

To  implement  the  asymptotic  analysis,  we  introduce  a  sequence  of  systems  S'„  indexed  by 
the  parameter  n  {n  —^  oo).  Each  system  has  the  same  structure  as  the  original  system  S. 
We  fix  {Cj}j>i?  ^  =  ^^---if^i  which  represent  the  interdemand  intervals,  and  allow  {^"'}j>i 
{?/"'} j>i  to  depend  on  n.  Thus,  the  parameters  associated  with  the  demand  c?,  and  i/^,-  do 
not  depend  on  n,  while  the  parameters  /i",  u^,  and  the  function  Fp{-),  i  =  1, . . . ,  A'  vary  with 


n.  To  emphasize  this  dependence  we  attach  an  upper  index  n  to  the  appropriate  processes 
yl",  T"  etc.  The  initial  inventories  C",  i  =  1, . . . ,  A'  with  respect  to  which  the  optimization 
procedure  will  be  implemented  can  also  depend  on  n. 

Let  p"  =  dil ji^  and  p"  =  ^j-xP".    Our  main  assumption  on  the  sequence  Sn  are  the 
following: 

^r^/i.,       K^^~^a^,    I   =    1 ,  •  •  •  , /T,  (3.1) 

riT  =  «"oS;  ,  (3.2) 

/j"  <  1  and    v/^(l  -  p")  =  1  .  (3.3) 

Observe  that  the  distribution  function  F"(<)  of  r/"'  is  given  by  P(7;"'  <  0  =  Fi{t/n)  = 
Ff{not/n),  where  F,(-)  is  defined  in  (2.7).  Assumption  (3.3)  is  the  heavy  traffic  condition 
for  the  sequence  Sn,  n  >  1,  which  requires  the  traffic  intensity  p^  to  converge  to  1.  As  a 
consequence  of  (3.1)-(3.3)  we  have 

d,  ^ 

p"^  —^  Pi  =  ^    as    n  — >  oo,  and    ^  p,  =  1  .  (3.4) 

/^'  .=1 

Assumption  (3.2)  implies  that  the  delays  are  increasing  as  n  — >  oo. 

Let  the  processes  Q^{-),  Z^{-),X^{-)  associated  with  system  Sn  be  defined  similarly  to 

(2.1)-(2.3): 

and 

Xnt)  =  Cr  -  Q^t)  -  Znt)  .  (3.5) 

For  the  processes  X",  Q",  Z"  and  the  cost  functional  associated  with  the  n^^  system,  we 
will  consider  two  types  of  rescaling:  Fluid  rescaling 

X"(0  =  -J^"(ni)  ,  Q^{t)  =  -Q'^int)  ,  Z"(0  =  -Z'^int)  ,  (3.6) 

n  n  n 

J(Ci",...,Q)  =  n~'limsup;J;£;  /    c{X''{t))dt  (3.7) 

T— oo      T       Jo 


I„  ,,e  remainder  of  tMs  sect.on,  we  descnbe     ^^    ^^^^  ^^^  ^^^^^^^^^  ^^^^^^  ^^^  ^^„,, 
processes  Q,  associated  wHh  product  ,. .  =  ^  ^^^^  ^^^^^  ,^„,  ,,,„,.  for 

;  ..e  ™acH.e  on  a  FCFS  U..  ";-    ^/^^^^^  ^^^  ^_„,.  «n,e  du.n,  ICl  .- 
n,.U-,dasssysten,s.Let   ^U  ,,fined  by 


n,.,Uiclass  systems,  .e.      w  defined  by 

the  production  mach,ne  ,s  .die.  Then        P 

h-  e's  queue  at  time  t.  Denote  the 
,e,ese„ts  the  total  atnount  of  «o..  P^;;;;;;;^^;;;:^  O^  =  "-'H-^'"''  '  *""■ 
.escaied  version  of  W^  by  ^"(0  =  "       _^^  ^^  ,,^  _,„ensionaf  teflect.on 

Fo.  any  funcfon  x(.)  Iron,  the  space  DlO.oo) 
capping  (see  Harrison  (1985));  that  is, 

Define  the  Brownian  motion  process  ^^^^^ 

WW  =  -'  +  «'.^w  +  ^' ''' '  '  - "  ' 

„«  with  variances  Sf^i-'kC'  "'"  + 
,       W^  and  W/  are  two  independent  Wiener  processes 
„here  W.   and  respectively.  Set 

Ef-,..^,<i^*i;'''3' »"'*''■''•     *  (3.11) 

B(.)  =  *(IV)()  , 

SO  that  B  IS  a  retlectea  d  ^^^^^ 

t=l 

.f  Theorem  1  of  Peterson  (1991). 
The  following  result  is  a  consequence  of  Theore 

8 


Proposition  3.1.     Under  assumptions  (3.1)-(3.3)  we  have 
where  — >  stands  for  weak  convergence  in  D^[0,oc). 

4      The  Limiting  Optimization  Problem 

As  mentioned  in  the  Introduction,  a  rather  crude  fluid  approximation  to  the  control  problem 
is  required  as  a  precursor  to  the  heavy  traflBc  diffusion-Gaussian  approximation. 
4.1.  Fluid  approximation.  For  brevity  of  notation,  we  introduce  G,(<)  =  1  —  Fi{t). 
Proposition  4.1.      Assume  (3.1)-(3.3)  and  n-^C^  ->  c\^\  i  =  1,...,A'.     Then  for  any 
C  >0, 


lim  sup  IX"  -  Xi(t)\  =  0    P  -  a.s.  , 


where 


Xi{t)  =  c^"^  -  d,  f  G,{u)ds  .  (4.1) 

JO 

Proof.  Proposition  3.1  implies  that 

Or(0^0,    n-^oo  ,  i  =  l,...,K  (4.2) 

uniformly  on  finite  time  intervals.  By  the  elementary  renewal  theorem  (Karlin  and  Taylor 
1975,  page  188),  n-'^D,{nt)  -^  d,t  as  n  ^  oo.  Since  0^(0  =  n-^[D,{nt)  -  A^{Tp{nt))],  we 
have 

n-'^A'^{Tp{nt))  ^  d,t  ,    n  ^  oo.  (4.3) 

Equation  (3.5)  impHes 

A^(Tr(nt)) 

xnt)  =  n-'cr  -  Qnt)  -  n-^   y:   w + ^r  >  "O  •  (4.4) 

i=l 

Observe  that  n~^Tf\  j  >  I  are  the  jump  moments  of  y4"(r"(n-))  and  n"^?;"*  =  tiq  Sj  bas 
the  distribution  function  F-'{nt)  =  F,(<).  Therefore,  the  last  term  in  (4.4)  is  the  rescaled 
number  of  customers  in  an  infinite  server  queue  with  service  time  distribution  i^(-)  and 


arrival  stream  A'^ {T-^ (nt)) .    Limit  theorems  for  such  models  were  considered  by  Borovkov 
(1984).  Theorem  1  in  Chapter  2  of  Borovkov  implies  that 


A?(Tr(nt))  t  t 

''       y       JivT  +  T?'  >  nt)  -^  d,  /   G,{t  -  u)du  =  d,  /   G,{u)du    P 

^  ■>  J  Jo  JQ 


uniformly  on  finite  time  intervals,  which  together  with  (4.2)  completes  the  proof.  □ 
Observe  that  the  cost  with  fluid  rescaling  can  be  written  as 

1  „   /-^ 

T— oo' 

Proposition  4.1  suggests  that  the  following  deterministic  optimization  problem  should  be 
considered  as  the  limiting  fluid  problem:  Choose  Cj    ,...  ,Cj^-   to 

1  „   fT 


J{C^, . . . ,  Q)  =  lim  sup  ^E  r  H{X^{t))dt 


minimize 

T-.00 


1       yj 
lize  iimsup— £^  /     c{xi{i))dt  ,  (4.5) 

T-.00    1      Jo 

where  a;,(-)  are  given  by  (4.1).  If  we  set 

cp)  =  d,  r  G,{u)du  ,  (4.6) 

then  lim(_oo  3:,(<)  =  0  and  the  associated  cost  is  equal  to  zero  and  is  minimal.  Since  c-     is 
the  limit  of  n'^C^  as  n  — >  cx),  we  set  C"  =  nc,-  '  +  AC",  where  AC"  =  o[n)^  and  turn  our 
attention  to  the  next  level  of  approximation. 
4.2.  Diffusion-Gaussian  approximation.    For  each  fixed  n,  define  the  processes  y"(<)  = 

y;"(0  =  ^."(0  -  K"(0  ,  <  >  0 ,  ^  =  l,...,/^,  (4.7) 

where 

V;"(<)  =  nc^P  -ndi  f   "  G,{u)du 
Jo 

and  c,-  '  are  given  by  (4.6).  Since  V-"(<)  — *  0  as  <  ^  cx)  for  each  n,  we  have 
i"(Cr,...,C^)  =  n-^/Mimsup^£;  /    //(X"(0)d<  =  n-^/Mimsup^£;  /    H{Y''{t))dt  . 

7'_oo      1        Jo  T-.00      1         Jo 

(4.8) 
10 


Thus,  our  goal  is  to  minimize  the  right-hand  side  of  (4.8). 
Consider  the  rescaled  process  ¥"{1)  =  n-^/^Y"-{nt),  t  >  0. 
Proposition  4.2  Assume  that  n~^/^AC"  — *  c,  ',  i  =  1, . . . ,  A'.    Then 


where 


Yi{t)  =  cp^  -  d,  f  Fi{t  -  u)dB{u)  -  /i.  /'  G,{t  -  u)dWf{u)  +  Oiit)  ,  (4.9) 

Jo  ^0 

and  Oi{-)  is  a  Gaussian  process  with  zero  mean  and  covariance  function 

E0,{t)e,{s)  =  d,  r  Fi{s  -  u)G,{t  -  u)du  ,    s  <t  (4.10) 

Jo 

independent  of  B{-),  Wf{-)  and  Oj{-),  j  /  i  {ij  =  1,. . . ,  A')- 

Proof.  Definitions  (3.5)  and  (4.7)  imply  that 

y;."(<)  =  n-v2Acr  -  n-^/2   Yl  A'/r'  +  ^r'  >  "<)  -  ^r(o +d,   G,it  -  u)du . 

Since  n~^T'^\  j  =  1,2, .. .  are  sequential  jump  moments  of  the  process  A^{TJ^{n-)),  for  any 
RCLL  function  /(•)  one  can  write 

K{Tr(nt)) 


E       fin-'r-)  =  f  f{u)dAnTnnu))  , 
j=i  •'^ 


which  is  a  Lebesgue-Stieltjes  integral  with  respect  to  a  nondecreasing  function.   Using  this 
fact  and  (3.2)  we  derive 

y;"(0  =  n-'/'AC."  -  f  F,{t-u)dQ'i[u)-  fG,{t-u)db^iu)  +  ent),        (4.11) 

Jo  Jo 

where 

A?iTp{nt)) 

0?{t)  =  n-'/'       j]      [/(noS]  +  n-V;-  <  t)  -  F,{t  -  n-'r^')]  .  (4.12) 

«=i 
The  process  6"  can  be  represented  as  an  integral  with  respect  to  a  random  field,  which  in  this 

case  is  a  stochastic  process  with  a  two-dimensional  variable.  To  derive  such  a  representation, 

we  define  the  random  fields 

[nt] 

UJ'it.x)  =  n"'/^  II[/(^.(«oSi)  <  a:)  -  max(x,  1)]  ,  ^  >  0  ,    x  >  0 

3=1 
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and 

Vrit,x)  =  n-'/'       Y.      [/(no^*  <  x)  -  F,(x)]  ,  i  >  0  ,    x>0. 

Hence, 

Vnt^x)  =  C/r(n-Mr(7;"(nO),F.(x))  .  (4.13) 

Observe  that  the  random  variable  F,(no^7/p  has  a  uniform  distribution  on  [0,1].  In  the 
statistical  literature  (see  Gaenssler  and  Stute  (1979),  for  example),  the  process  t/"  is  re- 
ferred to  as  the  stochastic  empirical  process.  The  trajectories  of  [/"(•,  •)  belong  to  the  space 
D^[0,oo)^  (see  Pardjanadze  and  Hmaladze  (1986)  and  Gaenssler  and  Stute),  which  is  the 
space  of  functions  u{t^x)  defined  on  [0,oo)  x  [0,oo)  such  that  the  limits 

u{t+,x+)  =    lim   u{s,y)  ,  u{t+,x-)  =    lim   u{s,y)  , 

slt,ylx  ait,y]x 

u{t—,x+)=    hm   u{s,y)  ,  u{t  —  ,x—)=    lim   u{s,y) 

s]t,ylx  s]t,y]x 

exist  and  u{t+,x+)  =  u{t,x).  The  space  i)^[0,oo)^  endowed  with  topology  similar  to  Sko- 
rohod's  Ji  topology  for  Z)[0,oo)  is  a  complete  metric  space  (see  Pardjanadze  and  Hmaladze 
for  details).  An  increment  of  a  function  u  €  D^[0,oo)^  is  defined  by 

□u((r,y),(/,x))  =  u{T,y)  -  u{t,x)  -  u{t,y)  +  u{t,x)  . 

Using  this  concept  one  can  define  the  space  of  functions  with  bounded  variation  in  the 
usual  way  (see  Kolmogorov  and  Fomin  (1975)).  Observe  that  the  trajectories  of  the  process 
(V^"(s,x),5  >  0,x  >  0)  belong  to  Z)^[0,oo)^  and  have  bounded  variation.  Therefore  the 
Lebesgue-Stieltjes  integral  below  is  well  defined  and  we  can  write 

e^{t)  =   f'  I{s  +  x<  t)dVi^{s,  t)  .  (4.14) 

Jo 

Using  the  definition  of  increments  given  above,  one  can  easily  check  that  (4.14)  is  equivalent 

to  (4.12).   Bickel  and  Wichura  (1971)  proved  that  the  empirical  process  Up{s,x)  converges 

weakly  in  D^[0,  oo)^  as  n  -^  oo  to  the  Gaussian  random  field  U{s,x)  with  zero  mean  and 

covariance  function 

EU{s,x)U{t,y)  =  sAt{x  Ay  -  xy)  ,    x  <  1,  j/  <  1  .  (4.15) 
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For  fixed  x,  U{-^x)  is  a  standard  Brownian  Motion,  and  U{s,-)  is  a  Brownian  Bridge  for 
fixed  s.  Since  t/"(-,)  is  independent  of  (5"()  and  £>"(•)  for  each  i  =  1,..,A',  Proposition 
3.1,    (4.3)  and  (4.13)  give  the  joint  convergence 

(Or(-),  A"(-),K"(-,-))  ^  {d,B{-),fl,Wf{-),V,i;.))  ,  n  -.  oo  ,  (4.16) 

where  Vi{t,x)  =  U{dit,Fi{x)).  To  complete  the  proof  we  need  to  take  the  Hmit  as  n  — +  oo 
in  (4.11).  Applying  the  integration  by  parts  formula  and  the  continuous  mapping  theorem, 
one  can  easily  verify  that  the  second  and  third  terms  in  (4.11)  converges  to  those  in  (4.9). 
Define 

Oi{t)  =  f  I{s  +  x  <  t)dV,{s,x)  , 

.70 

where  the  integral  is  interpreted  as  the  mean-square  integral  with  respect  to  a  Gaussian  field. 
Then  (4.14)  and  (4.16)  imply  that 

e^  ^e,  ,    n^oo  .  (4.17) 

Unfortunately,  the  integration  by  parts  formula  cannot  be  applied  directly  in  this  case,  since 
I{s  +  X  <  t)  as  a.  function  of  {s,x)  does  not  belong  to  D^[0,oo)^.  Hence,  a  careful  analysis 
is  required  to  justify  this  convergence.  We  do  not  provide  here  the  proof  of  this  result,  since 
it  was  carried  out  in  Theorem  2  of  Borovkov  by  different  techniques.  Readers  are  referred  to 
Krichagina  and  Puhalskii  (1993)  for  an  alternative  proof  based  on  the  representation  (4.14). 
The  joint  convergence  of  (4.16)  and  (4.17)  allows  us  to  take  the  hmit  as  n  — >  oo  in  (4.11)  and 
obtain  the  desired  result.  It  follows  from  (4.15)  that  the  covariance  function  Ki{u,v,x,y)  of 
the  process  K(-,  •)  =  Ui{dit,  Fi{x))  is  given  by 

Ki{u,v,x,y)  =  EVi{u,x)Vi{v,y)  =  di{u  A  v){F{x)  A  F(y)  -  F{x)F{y))  . 

The  properties  of  the  mean-square  integral  imply  that 

Jroo     TOO  ^oo 

'      /     I{u  +  x<s)dV,{u,x)         Iiv  +  y<t)dV,{v,y) 
0      ^0  Jo 

roo     roo     ^oo     roo 

=     /      /      /      /     I{u  +  x<s)I{v  +  y<t)dKi{u,v,x,y). 
Jo     Jo     Jo     Jo 
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Since  duAv  =  I{u  =  v)du,  dF,{x)AF,{y)  =  I{x  =  y)dF,{x)  and  dFi{x)F,iy)  =  dF,{x)dF,{y), 
the  right-hand  side  ior  s  <  t  is  equal  to 

/■OO        /■OO  /-co        /•TO        /-OO 

(fj/      /     I{u  +  x<s)dudF,{x)-  /      /     /(ii  +  a-<5)/(w  +  j/<0<iuc?F,(x)(fF.(r/)}  . 

Jo     Jo  Jo     Jo     Jo 

Calculating  the  integrals  we  obtain  (4.10).  □ 

For  each  fixed  i,  let  F,{t),  i  <  0  be  the  distribution  on  (  — oo,0]  generated  by  F,(-)  as 

follows: 

F,{t)  =  1  -  F{-t)  =  G,{-t)  ,    t<0  . 

Denote  by  'yf  and  'yf^'  (the  superscripts  represent  demand  and  delay,  respectively)  two  in- 
dependent Gaussian  random  variables  with  zero  means  and  variances  diV^-  f^  G^{u)du  and 
di  f^  Fi{u)Gi{u)du,  respectively.  Let  5*'(-)  be  the  stationary  RBM  with  the  same  structure 
and  parameters  as  B{-)  (see  (3.10)-(3.12)).  Since  B^*  is  a  stationary  process  we  can  assume 
that  it  is  defined  for  /  ^  (— oo,+oo).  Define 

7f  =  d,  f    B'\u)dF{u)  ,  (4.18) 

J  —  OO 

which  is  the  stationary  Brownian  motion  averaged  by  the  distribution  function  Fi{t)   = 
Gi{—t).  The  following  proposition  describes  the  steady-state  distribution  of  yi(-). 
Proposition  4.3.     Under  the  assumptions  and  definitions  above,  Yi{t)  — >  5^(00)  as  t  —^  00, 
where 

y;(oo)  =  cp^-7f-7f  +  7-^'.  (4.19) 

Proof.  To  prove  the  result  we  take  the  limit  in  (4.9)  as  <  — >  00.  Our  first  step  is  to  show 
that 

R,{t)  =  r  F,{t  -  u)dB{u)  4  /     B''{u)dF,{u)  .  (4.20) 

Jo  J- 00 

An  application  of  the  integration  by  parts  formula  yields 

R-{t)  =  /'  FAt  -  u)dB{u)  =  -  f  B{u)dFi{t  -u)=  f  B{u)dG,{t  -  u)  .  (4.21) 

Jo  Jo  Jo 

For  yV  such  that  0  <  TV  <  <,  define 

R'^{t)  =  f     B{u)dG,it  -  u)  .  (4.22) 

Ji-N 
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By  Theorem  4.2  of  Billingsley,  (4.20)  will  follow  from  the  statements  below: 

R^{t)  -^  I     B''(u)dF,{u)  ,  t^+oo   for  fixed  N  ;  (4.23) 

J-N 


L 


°   B'\u)dFi{u)  -i  f    B''\u)dF,{u)  ,    AT  -^  +00  ;  (4.24) 

J  —  oo 


-N 

N 


lim  \im sup P{\R,{t)-R^{t)\  >  e)  =  0  ,    e>0.  (4.25) 

N^co      t— ►oo 

To  prove  (4.23),  we  start  by  deriving  R^{t)  =  J^  B{t  -  N  +  u)dG,{N  -  u)  from  (4.22). 
Standard  arguments  can  be  used  to  show  that  B{t  +  ■)—>■  B^\-)  as  <  — >  cx).  For  fixed  A^, 
the  continuous  mapping  theorem  gives  ^f  (<)  — >  /q  B"\u)dGi{N  —  u),  t  —^  oo.  Since  B'* 
is  a  stationary  process,  we  have  /q  B'\u)dGi{N  —  u)  =  J_^  B'\u)dGi{—u),  and  (4.23)  is 
proved. 

Because  EB'\u)  =  const,  it  follows  that  E  jZ^  B'\u)dF,{u)  =  {const)Fi{-N)  = 
{const)GiiN)  ->  0  as  TV  ^  oo.  Thus,  (4.24)  is  verified. 

Equations  (4.21)  and  (4.22)  imply  that  R,{t)-R^{t)  =  f^'^  B{u)dGi{t-u).  The  closed 
form  expression  for  the  distribution  of  B{t),  t  >  0  (see  p.  15  of  Harrison  1985)  implies  that 
EB{t)  <  const.  Thus  E\R,{t)  -  R^{t)\  <  {const){G,iN)  -  Gi{t)).  Applying  Chebyshev's 
inequahty  and  allowing  i  — >  oo  and  then  N  —^  oo  yields  (4.25),  and  hence  (4.20). 

For  each  i  >  0,  the  random  variables  //;  Jq  Gi{t  —  u)dWf{u)  and  0i{t)  are  Gaussian  with 
zero  mean  and  variances  dii^j-  jQG]{u)du  and  (f, /q  F,(u)G,(u)c?u,  respectively.  As  <  — >  oo, 
they  converge  to  ff  and  7^^',  respectively.  The  joint  convergence  follows  from  the  structure 
of  B  and  the  independence  of  $i  and  B,  Wf.  Thus,  one  can  take  the  limit  in  (4.9)  as  i  ^  00 
and  obtain  the  result.  □ 

The  cost  functional  given  by  (4.8)  can  be  expressed  in  terms  of  ¥"■  as 


limsup;^;^  /    H{Y''{s))ds  . 
r-.oo    J      Jo 


Based  on  this  fact  and  Propositions  4.2-4.3,  the  limiting  optimization  problem  corresponding 
to  the  diffusion-Gaussian  approximation  is  to  choose  c[   , . . . ,  c^-   to 


mmimize  lim 

T— »oo 


sup  ^E  [    H{Y{s))ds  =  EH{Y {00))  (4.26) 

.00    -I      Jo 
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where  V'(oo)  =  (y'i(oo), . . . ,  ^^-(oo))  is  given  by  (4.19).  The  additive  structure  of  the  cost 
and  steady-state  equalities  (4.19)  imply  that  the  minimization  can  be  implemented  for  each 
component  i  individually.  Thus,  (4.26)  is  equivalent  to  finding 

min  EH,{Y,{oo))  (4.27) 

for  each  i  =  1, . . . ,  K .  Let  c|  ,  z  =  1, . . . ,  A'  be  the  solution  to  (4.27).  Then  our  two-stage 
asymptotic  analysis  implies  that  nearly  optimal  values  for  the  parameters  C",...,C^-  are 
given  by 

C^  =  ncl'^  +  V^c*''  ,  i  =  l,...,K  .  (4.28) 

5      Two-Stage  Approximation  for  the  Original  System 

In  this  section,  we  employ  the  limiting  fluid  and  diffusion-Gaussian  problems  as  approxima- 
tions for  the  original  optimization  problem,  and  provide  numerical  solutions  under  certain 
additional  assumptions.  The  relation  between  the  original  system  S  and  the  sequence  Sn, 
n  >  1  is  that  the  dynamics  of  S  can  be  identified  with  that  of  S'no,  where  no  is  given  by  (2.6), 
if  we  set  //,  —  ^i^°  and  i/^,  —  {i^aiY-  Notice,  however,  that  the  hmiting  parameters  /i,  and 
j/^-,  i  —  1,. . .  ,K  are  not  directly  connected  with  those  of  the  original  system  S;  they  were 
introduced  to  justify  the  Hmiting  procedure  in  Propositions  3.1  and  4.1-4.2.  Therefore,  to 
obtain  the  diffusion-Gaussian  approximation  for  the  original  system  5,  we  substitute  limiting 
"tilde"  parameters  by  the  corresponding  parameters  of  the  original  system. 

The  original  optimization  problem  with  the  cost  functional  (2.4)  normalized  by  the  factor 
Uq^  can  be  approximated  by  the  fluid  model  (4.1),  (4.5).  As  in  (4.6),  the  optimal  parameters 
are 

cS'^  =  d,  r  G,{u)du  =  T^  ,  z  =  1, . . . ,  A'  . 
Jo  A,  no 

In  accordance  with  (4.28),  the  nearly  optimal  parameters  at  the  fluid  level  of  approximation 
are 

C.«nocS')  =  ^.  (5.1) 
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Since  d  are  the  initial  inventories,  this  expression  admits  the  following  explanation.  The 
system  operates  in  a  balanced  regime  and  the  average  production  rate  is  equal  to  the  average 
demand  rate  d,  for  each  product  i.  Therefore,  the  number  of  units  in  the  delay  node  plus  the 
number  in  inventory  is  approximately  equal  to  the  initial  inventory  C,  for  any  i  =  1,. . .  ^K. 
Over  the  long  run,  the  average  number  of  product  i  units  at  the  delay  node  is  equal  to  di/\i, 
and  therefore  product's  Vs  inventory  is  roughly  C,  —  di/\,.  If  C,  differs  from  (5.1),  then 
there  would  be  a  constant  backlog  or  a  constant  excess  of  inventory  over  the  long  run,  which 
would  increase  the  cost  functional. 

The  diffusion-Gaussian  approximation  corresponds  to  the  original  problem  with  the  cost 
functional  (2.4)  normalized  by  the  factor  no  .  To  find  the  corresponding  term  y/noc-  in 
(4.28),  we  must  solve  the  minimization  problem 

min  EH,{Yi{oo))ds  (5.2) 

for  each  i  =  !,...,/<',  where 

K(cx))  =  c<^)-7f-7f  +  7f^  (5.3) 

7f  -  d,  f    B'\u)dF,{u)  , 

J  —  CO 

B^^  is  the  stationary  RBM  on  [0,  oo)  with  drift  —p  and  variance  X^aLi  ^ArA'fc  ^('^ait  +  '^dk)i 
and  ff  and  ff^'  are  two  independent  Gaussian  random  variables  with  zero  mean  and  vari- 
ances diU^-  f^  G^{u)du  and  di  J^  F,{u)Gi{u)du,  respectively.  Once  the  parameter  c-  '  that 
minimizes  the  cost  (5.2)  is  found,  we  use  (4.28)  to  set 

C,  =  ci'^no  +  cj"v/n^  =  d,/X,  +  c!'V(l  -  p)  ■  (5-4) 

The  fluid  approximation  yields  the  crude  deterministic  description  of  the  system  in  terms 
of  the  average  values,  and  the  diffusion- Gaussian  approximation  provides  the  distribution  of 
fluctuations  around  these  average  values.  Thus,  the  second  term  on  the  right-hand  side  of 
(5.3)  characterizes  the  steady-state  inventory's  fluctuations  caused  by  the  variation  in  the 
size  of  the  machine's  queue.  The  third  term  represents  the  impact  of  demand  variabihty,  and 
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the  fourth  term  quantifies  the  effect  of  fluctuations  in  delay.  By  (2.7),  F,(-)  is  the  resettled 
distribution  function  of  the  delay  times,  and  it  plays  an  integral  role  in  the  definition  of 
these  quantities.  Due  to  the  rescaHng  procedure  {t  — >  nt),  we  observe  all  processes,  as  well 
as  the  length  of  delays,  in  the  compressed  time  scale.  Thus,  our  analysis  is  oriented  towards 
systems  with  long  and  significantly  uncertain  delays  (that  is,  "wide"  distribution  functions 

Unfortunately,  problem  (5.2)  does  not  admit  an  exphcit  solution.  The  main  difficulty 
is  that  the  random  variables  7f  and  ^f  are  dependent  and  their  joint  distribution  is  not 
easily  derived.  However,  the  expression  for  the  steady-state  distribution  (5.3)  reveals  the 
structure  of  the  limiting  behavior  of  production-distribution  systems  with  long  delays,  and 
can  be  used  for  their  analysis.  To  illustrate  the  use  of  our  analysis,  we  assume  for  simplicity 
that  F^  is  a  uniform  distribution  on  the  interval  [(1  -  af''')/\,,  (1  -f  af^')/A,],  where  af^'  <  1 
characterizes  the  uncertainty  of  delays.  In  this  case,  the  rescaled  distribution  function  Fi  is 
uniform  on  an  interval  of  length  /,  =  2af^'{l  —  p)^/\,,  and  ff  is  equal  in  distribution  to 

^  /'■  B'\u)du  =   C  d,B'\hu)du  .  (5.5) 

/;    Jo  Jo 

Let  B-*{-)  be  the  unitless  RBM  on  [0,  oo)  with  drift  —p  and  variance 

Define  the  unitless  constant  I,  =  d,h  =  2af''{l  -  pYidJX,).    Since  d,B'\l,-)  =  Bf{Li-), 

equation  (5.5)  gives 

7f  =  /'  Bf{L,u)du  =  ]-  r  Bf{u)du  .  (5.6) 

Jo  Li  Jo 

The  variances  of  ff  and  jf^'  are  equal  to 

4L,/3   and     L./6,  (5.7) 

respectively.  In  the  remainder  of  this  paper,  we  investigate  three  cases  according  to  the 
magnitude  of  Z/,. 

Case  1.  Suppose  that  the  parameter  L,  is  small.  This  happens  when  the  delays  are  not  long 
enough  or  have  low  variability,  or  the  demand  rate  is  small.  Then  the  distribution  of  7,-    is 
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close  to  that  of  Bf'(O),  which  is  exponential  with  parameter  'Ip/E^.    By  (5.7),  the  random 
variables  7/  and  7^^'  have  sniall  variances  and  can  be  neglected.  It  is  easy  to  calculate  that 


(2)       ■^^  1     /  h,  +  b, 


V;2 


ip     -       h, 

minimizes  EH,{c\''  —  jf).    Thus,  for  small  L,,  the  variability  in  the  size  of  the  machine  s 
queue  has  the  largest  impact  on  the  fluctuations  of  the  inventory. 
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Figure  2.  The  simulated  distributions  of  f^  for  various  value.s  of  L,. 

Case  2.   Suppose  that  L,  is  large.  By  (5.7),  the  variances  of  7/  and  -ff^'  are  large  a?  well. 
By  the  Birkhoff- ivhinchinc  theorem,  we  have 

--  /""  Bfiu)du  ^  EBf  {00)  =  ^  ,    I,  ^  00  . 
L,  Jo  Ip 

Hence,  by  (5.6),  the  distribution  of  ■^f  for  large  //,  has  small  variability,  and  one  .,an  let 

7,^  Ri  TJlj2p.  Since  '^f  and  7,*^^'  are  independent,  their  sum  has  a  Gaussian  distribati:'n  with 
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zero  mean  and  variance  Li{i/jj3  +  1/6).  The  value  of  the  parameter 

cS''  =  argmin,//.(c-S  +  ^f  +  ^_<^e/) 

2p 

is  easy  to  find  numerically. 

Case  3.  Problem  (5.2)  is  more  difficult  when  L,  is  neither  small  nor  large.  To  make 
further  progress  in  this  case,  we  make  the  additional  assumption  that  'yf  «  0  and  can  be 
neglected.  That  is,  i/J,  is  small  and  the  demand  is  nearly  constant.  As  alluded  to  earlier,  this 
assumption  is  required  because  jf  and  7,^  are  dependent  random  variables.  The  distribution 
of  'yf  also  cannot  be  calculated  explicitly.  However,  a  one-dimensional  reflected  diffusion  is 
easily  simulated. 

The  simulated  distributions  of  jf  for  various  values  of  Li  are  presented  in  Figure  2.  As 
predicted  by  Case  1,  we  observe  that  for  small  values  of  L,  the  distribution  is  close  to  the 
exponential  with  parameter  2p/S^  shown  by  the  dashed  line.  In  accordance  with  Case  2, 
Figure  2  also  indicates  that  for  larger  values  of  L,  the  distribution  is  concentrated  around 
the  mean  of  the  exponential  distribution  shown  by  the  vertical  line.  By  calculating  the  cost 
function  Hi{x)  averaged  by  the  convolution  of  the  distribution  of  7P  found  by  simulation 

('2\ 

and  a  normal  distribution  with  mean  c\  and  variance  Li/G,  we  were  able  to  obtain  the  cost 
for  any  given  c,  and  find  the  minimum.  Tables  1  and  2  compare  our  results  (COST  TH) 
with  the  results  obtained  via  simulation  of  the  original  system  (COST  SIM)  for  one-  and 
two-dimensional  problems.  The  service  times  were  assumed  to  be  uniformly  distributed  on 
[(1  —  Q,)///i,  (1  -f  Qi)///,],  i  =  1, . . . ,  A'.  In  all  four  problems,  our  cost  minimizing  values  of 
Ci  are  very  close  to  those  obtained  via  simulation. 
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Problem  1  parameters:  Hi  = 
Number  of  simulated  jobs  : 

=  0.9-* 

2.5  X 

1.0,ai 

=  0.8, 

a^^'  = 

0.8,1/Ai  =50.0. 

Ci 

50.0 

51.0 

52.0 

5.3.0 

54.0 

55.0 

56.0 

57.0 

58.0 

59.0 

COST  SIM 

10.41 

8.. 36 

6.94 

6.09 

5.77 

5.89 

6.32 

7.00 

7.80 

8.70 

COST  TH 

11.15 

8.79 

7.09 

6.03 

5.52 

5.48 

5.80 

6.38 

7.15 

8.03 

Problem  2  parameters:  c?i  =  1.0, 1/A 
'*^umber  of  simulated  jobs:  5.0  x  10^ 

=  100.0,/ii  =0.9-\a 

=  0.8,0"^"'  =  0.5. 

Ci 

100.0 

101.0 

102.0 

103.0 

104.0 

105.0 

106.0 

107.0 

108.0 

109.0 

110.0 

COST  SIM 

11.29 

9.30 

7.84 

6.96 

6.57 

6.59 

6.93 

7.50 

8.23 

9.07 

9.88 

COST  TH 

12.09 

9.76 

8.04 

6.88 

6.24 

6.04 

6.22 

6.67 

7.33 

8.13 

9.03 

Problem  3  parameters:  di  =  1.0, 1/A 
Number  of  simulated  jobs:  7.5  x  10^ 

1  =  100.0,^1  = 

0.9-\o 

.  =  0.8, 

a''^'  =  0.99. 

Cr 

100.0 

101.0 

102.0 

103.0 

104.0 

105.0 

106.0 

107.0 

108.0 

109.0 

110.0 

COST  SIM 

14.94 

12.97 

11.41 

10.25 

9.46 

9.01 

8.86 

8.98 

9.31 

9.82 

10.47 

COST  TH 

15.8 

13.59 

11.79 

10.40 

9.41 

8.78 

8.48 

8.46 

8.69 

9.14 

9.70 

Table  1:  Three  one-dimensional  problems. 
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Parameters:  /Xi  =  1.0,  //2  =  0.8,  d^  =  0.5,  (ij  =  0.32,  {p  =  0.9),  l/A^  =  I/A2  =  50.0, 

Qi  =0.6,a2  =  0.8,  af  = 
Number  of  simulated  jobs:  2.0  x  10^. 


Q^ '  =  0.8. 


COST  SIM: 

C2         Cl 

21.0 

26.0 

27.0 

28.0 

29.0 

30.0 

35.0 

11.0 

53.20 

34.36 

32.85 

32.17 

32.17 

32.64 

37.24 

16.0 

31.47 

11.00 

11.12 

10.44 

10.44 

10.90 

15.49 

17.0 

29.13 

10.47 

8.96 

8.28 

8.28 

8.76 

13.36 

18.0 

28.27 

9.43 

7.91 

7.24 

7.24 

7.71 

12.31 

19.0 

28.14 

9.30 

7.78 

7.11 

7.11 

7.58 

12.18 

20.0 

28.62 

9.78 

8.27 

7.59 

7.59 

8.06 

12.67 

25.0 

33.35 

14.51 

13.00 

12.32 

12. .32 

12.79 

17.40 

COST  TH: 


C2         Cl 

21.0 

26.0 

27.0 

28.0 

29.0 

30.0 

35.0 

11.0 

48.47 

31.31 

30.32 

30.00 

30.40 

31.06 

35.89 

16.0 

27.60 

10.44 

9.45 

9.21 

9.52 

10.19 

15.02 

17.0 

25.88 

8.72 

7.73 

7.50 

7.80 

8.74 

13.30 

18.0 

25.25 

8.09 

7.10 

6.86 

7.17 

7.84 

12.67 

19.0 

25.41 

8.25 

7.26 

7.02 

7.33 

8.00 

12.83 

20.0 

26.06 

8.90 

7.91 

7.60 

7.98 

8.65 

13.48 

25.0 

30.93 

13.77 

12.78 

12.60 

12.85 

13.52 

18.30 

Table  2:  A  two-dimensional  problem. 
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